/*-----------------------------------------------------------------------------\
 * |                         Matrices Header file
 * |
 * | This is header file for module matrices, which defines basic linear
 * | algebra objects and implements basic algebraic functionality.
 * |
 * | Dependencies: needs definition of double and LS_ACCURACY
 * \-----------------------------------------------------------------------------*/

#ifndef __MATRICES__
#define __MATRICES__ 1

#include "glodef.h" /* provides definition of double and LS_ACCURACY */

/*---------------------------------------------------------------------------\
 * | Definitions of MATRIX and VECTOR types and declarations of routines that
 * | provide the memory management for these types
 * \---------------------------------------------------------------------------*/

typedef double *VECTOR;
typedef double **MATRIX;

VECTOR  allocV(int dim);   /* allocate memory for vector */
void freeV(VECTOR *V);                 /* release VECTOR */
/* alloc memory for MATRIX */
MATRIX  allocM(int rows, int cols);
void freeM(MATRIX *M); /* release MATRIX */

void CopyVector(const double *S, VECTOR T, int Length); /* copies S into T*/
void CleanVector(VECTOR V, int length); /* sets entries in vector V to 0 */

/*---------------------------------------------------------------------------\
 * | Sparse_Matrix type and routines
 * \---------------------------------------------------------------------------*/

// #define USE_SPSOLVER
// #define USE_HYPRE

#if defined(USE_SPSOLVER)
/* you must use CSC with SPSOLVER */
#define USE_CSC 
#else
/* you may or may not use CSC here. not using CSC is marginally faster */
/* #define USE_CSC  */
#endif


typedef struct _SPARSE_MATRIX
{
    int rows;           /* Number of rows in matrix (first dimension) */
    int cols;           /* Number of columns */
    int Symmetric;      /* Indicates if the matrix is symmetric */
    char data_type[3];
    void *SPARSE_DATA;  /* Pointer to internal data structure that contains
     * the nonzero elements of the matrix and other
     * useful information. User subroutines defined below
     * can be used to access this data. For definition and
     * description of this structure see remarks in
     * matrices.c */
} Sparse_Matrix;

void PrintVector(double *V, int dim, char *filename);
void PrintSMatrix(Sparse_Matrix *M, char *filename);

/*  The following routines allocate and initialize Sparse_Matrix structures. */
/*  Returns TRUE if successful, or FALSE if error occurs */
int CreateSM(Sparse_Matrix *M, int rows, int cols, int sym, char *type);
/*  To release memory used up by a Sparse_Matrix use: */
void DestroySM(Sparse_Matrix *M);

void SM_BeginStructure(Sparse_Matrix *M);
void SM_EndStructure(Sparse_Matrix *M);
void SM_BeginAssemble(Sparse_Matrix *M);
void SM_EndAssemble(Sparse_Matrix *M);

/*  SM_SetEntry assigns the given value to an entry in the matrix */
void SM_SetEntry(Sparse_Matrix *M, int row, int col, double AValue);
/*  SM_AddToEntry incteases an entry in the matrix with a given value */
void SM_AddToEntry(Sparse_Matrix *M, int row, int col, double AValue);

/*  SM_GetEntry returns the value of an entry in the matrix */
double SM_GetEntry(Sparse_Matrix *M, int row, int col);

/* Sets all entries of matrix to 0. */
void CleanSparseMatrix(Sparse_Matrix *M);

/****************************************************************************
 *  The following routine "visits" every nonzero entry in the matrix and    *
 *  performs certain "action" on it. The "action" is a  user defined        *
 *  subroutine that will be envoked for each nonzero element and will be    *
 *  provided with four pieces of information. Namely, these are the row and *
 *  column coordinates of the entry, its value and a pointer "Data". The    *
 *  user can use this pointer to pass custom data to the "action"-routine.  *
 ****************************************************************************/
typedef int (*SM_ENTRY_ACTION)(int row, int col, double Val, void *Data);
int ForEachNonZeroEntryDo(Sparse_Matrix *M,
        SM_ENTRY_ACTION action, void *Data);

/*  The following routine quickly multiplies matrix by vector */
void MultiplySparseByVector(Sparse_Matrix *M, int dim_vect,
        const double *vect, int dim_res,  double *res);
void MultiplySparseTransByVector(Sparse_Matrix *M, int dim_vect,
        const double *vect, int dim_res,  double *res);
/*  This routine calculates a linear combination of n matrices. The matrices */
/*  are in the array M and the corresponding coefficients in array C. */
void AddSparseMatrices(int n, double C[], Sparse_Matrix *M[],
        Sparse_Matrix *Res);

/*-----------------------------------------------------------------------------\
 * | Linear System Solvers
 * \-----------------------------------------------------------------------------*/

/*  This constant defines the target accuracy of the iterative solvers */
#define SOLVER_TOLERANCE      LS_ACCURACY

/*  This will solve system with regular matrix */
int SolveLinearSystem(const MATRIX m, int dim, const double *rhs, VECTOR sol);
/*  This is an interface to routine that uses Conjugate Gradient Method */
/*  applied on a system with Sparse Matrix */
/*
  int SolveSparseSystem(Sparse_Matrix *M, int dim, VECTOR rhs, VECTOR *sol,
  double *err, double *res);
 */

/********** The conjugate gradient method ***********************************/

typedef void (*matr_prod_t)(const double *X, double *F, int trans, void *mp_data);
typedef void (*prec_solv_t)(double *X, const double *F, void *ps_data);

/* General PCG solver. */
void solve_cgm(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU);

void solve_bcgm(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU);

/* CG for sparse matrices.  If matrix is not symmetric uses Bi-CG method. */
void sparse_cgm(Sparse_Matrix *M, double *X, const double *F, double TOL);

/* Preconditioned CG for sparse matrices with diagonal preconditioner */
void sparse_pcgm_diag(Sparse_Matrix *M, double *P, double *X,
        const double *F, double TOL);

/* NOTE: preconditioning for the bi-cg method is not available at this time. */

/* general sparse solver */
void SolveSparse(Sparse_Matrix *M, double *X, const double *F, double TOL);

#endif

